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Abstract 

We are developing a framework for multiscale computation which 
enables models at a "microscopic" level of description, for example 
Lattice Boltzmann, Monte Carlo or Molecular Dynamics simulators, 
to perform modelling tasks at the "macroscopic" length scales of in- 
terest. The plan is to use the microscopic rules restricted to small 
patches of the domain, the "teeth" , followed by interpolation to esti- 
mate macroscopic fields in the "gaps". The challenge addressed here 
is to find general boundary conditions for the patches of microscopic 
simulators that appropriately connect the widely separated "teeth" to 
achieve high order accuracy over the macroscale. Here we start explor- 
ing the issues in the simplest case when the microscopic simulator is 
the quintessential example of a partial differential equation. For this 
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case analytic solutions provide comparisons. We argue that classic 
high-order interpolation provides patch boundary conditions which 
achieve arbitrarily high-order consistency in the gap-tooth scheme, 
and with care are numerically stable. The high-order consistency is 
demonstrated on a class of linear partial differential equations in two 
ways: firstly using the dynamical systems approach of holistic dis- 
cretisation; and secondly through the eigenvalues of selected numer- 
ical problems. When applied to patches of microscopic simulations 
these patch boundary conditions should achieve efficient macroscale 
simulation. 
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1 Introduction 



As a preliminary illustration of the gap-tooth scheme consider sim- 

ulating the diffusion and nonlinear advection of the viscous Burgers' equation 

du iqq u ® u — ® 2u m 
dt dx dx 2 
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Figure 1: Gap-tooth solution of Burgers' equation (JTJ) on [0, 2tc] through 
microsimulation on 8 teeth, each of small width 7r/20; the teeth are coupled 
by special patch boundary conditions. 

Suppose our aim is to simulate the evolution of fields u(x, t) periodic in x on 
the macroscopic length scale 2tt. See in Figure Q the continuous time evolu- 
tion on m = 8 grid "points" in space with macroscopic spacing H = it/A. 
However, each "point" is actually a microscopic patch of width h = 7r/20. 
Further, the only knowledge that the macroscopic evolution has of Burgers' 
pde (Q) is through the detailed simulation of the pde within each patch; 
here we obtain this local detailed simulation via a discretisation of (|T|) on a 
microscopic spatial grid of n — 11 points within each patch, Ax = 0.0175, 
and on a microscopic time step of At ~ 10 -4 . This fine scale discretisation 
of Burgers' pde (0) represents a finely detailed model or particle simulation 
that is too expensive to use over the entire macroscopic domain. Our task 
here is to begin to show how such microscopic simulations in relatively small 
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Figure 2: Gap-tooth solution of Burgers' equation (0) on [0, 2-n] on 8 teeth 
each of small width 7r/20 and coupled by special patch boundary conditions. 
Solutions u(x,t) + At are plotted at five times t = : 0.025 : 0.1 in different 
colours and connected by yellow lines and with the vertical displacement of At 
to help distinguish the plots. 



patches of space may be coupled by patch boundary conditions derived in 
Section |21 to ensure high order accuracy over the macroscopic domain. 

The example of Figures Q and 121 shows us that there are two time scales 
in the simulation. Rapidly, the initial internal structure within each tooth 
(black curves in Figure EJ) smooths by diffusion on the microscopic time- 
scale to a local quasi-equilibrium (blue curves). Then, over longer times the 
inter-patch coupling exchanges information between the teeth to guide how 
the local quasi-equilibria evolve over macroscopic times. See the dynamics 
in Figure |2J the broad hump initially centred around x ~ 4 is nonlinearly 
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advected to the right to wrap around to about x ~ 1 at the end (red curves); 
whereas the short hump initially at x ~ 1 is dissipated quickly against the 
slow moving region at x « 2 . Because of the two time scales, we plan future 
research to implement "coarse grained" integration 011] which uses just short 
bursts of microscopic integration to then extrapolate over a macroscopic time 
step. The result will then be a scheme where the microscopic simulations are 
only needed for relatively small patches in space-time. However, here we 
concentrate on only the issue of the macroscopic coupling of small patches 
across space. 

The method of "holistic discretisation" , developed by Roberts and Macken- 
zie [TUl HU C21 Hj j creates discretisations on a macroscopic grid using system- 
atically obtained analytic approximations for the subgrid field. The analytic 
solutions of Section El in this method are analogous to the microscopic system 
simulators in the gap-tooth scheme: they both provide microscopic solutions 
which are macroscopically coupled to neighbouring elements. This dynami- 
cal systems approach adapts the patch boundary conditions to support the 
modelling by centre manifold theory Ilj. Then the equivalent pde of the 
macroscopic dynamic model is found, in Section EJ to confirm the high order 
consistency for a wide class of linear pdes. 

Lastly, in Section 01 we consider a numerical time integrator for the dif- 
fusion equation on patches. The eigenvalues of the integrator again confirm 
the high order accuracy of the proposed patch boundary conditions. 

2 Couple the patches 

In this section we develop a coupling of the internal dynamics of patches with 
their neighbours to achieve high order consistency. We construct a boundary 
condition for the flux on the edge of the microscopic patches that is a natural 
interpolation of the surrounding macroscopic field. The number of required 
boundary conditions will depend upon the microscopic simulator |H1 e.g.], 
and possibilities other than the flux remain to be explored. But here we 
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know that boundary conditions on the flux should give rise to well-posed 
diffusion-like dynamics. 

We introduce the notation in which we typically use capital letters for 
macroscopic quantities and lower case letters for microscopic quantities. Thus 
let each of m patches be centred on equi-spaced grid points x = Xj = jH 
seen in Figures [0 and 121 Let each patch be of width h. Then the edge of 
a patch is a distance h/2 from its grid point, a fraction r = h/(2H) to the 
neighbouring grid point: when r = | the neighbouring patches meet and 
there would be no gap, as in holistic discretisation (THj. Here we expect 
the fraction r to be small so that the patches are a relatively small part of 
the physical domain. For example, r = 1/10 in Figures Q and 121 Now let 
Vj(x,t) be the microscopic field in the jth patch. 

We use the following identities for discrete operators [8. on a step size of 
the macroscopic grid and are careful whether we are using as a step of H 
in x or a step of 1 in j. In terms of the shift operator, Ev(x, t) = v (x + H, t) 



or equivalent ly EUj = Uj + i. 

centred mean fi = \(E 1/2 + E~ 1/2 ) , (2) 

centred difference 5 = E l/2 - E~ 1/2 , (3) 

translate/shift E = 1 + /j,5 + \5 2 , (4) 

derivative in x Hd x = 2 sinh" 1 \8 , (5) 

an identity ft 2 = 1 + \5 2 . (6) 



For example, the derivative of the microscopic field on the edge of a patch, 
at (Xj±rH, t), may be obtained from Vj through applying the operator 

E ±r Hd x = (l + n8+ \5 2 ) ±r 2 sinh" 1 \8 by © and © 

= [l±rfi5 + 0(5 2 )][5 + 0(5 3 )] 

= 5±rfi5 2 + 0(5 3 ) 

= fi5±r5 2 + (6 3 ) by©. (7) 

This last operator just involves evaluation at the grid points Xj and hence 
is evaluated from the macroscopic grid values Uj. This provides the same 
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approximation for the microscopic gradient as obtained by quadratic interpo- 
lation through the neighbouring macroscopic grid values |5J e.g.]. We proceed 
to modify such a patch boundary condition in order to obtain higher order 
consistency with the surrounding macroscopic variations. 

For arbitrary order consistency, as the macroscopic grid size H — > or 
as the gradients become small, repeat the previous analysis but retain more 
terms, and using (0) to replace /i 2 terms: 

E ±r Hd x = (1 + //<* + i5 2 ) ±r 2sinh _1 i5 

= l i5±r5 2 - (i - ±r 2 )/i<5 3 T - ±r 2 )5 4 

^30 8' ^24' Ir ^'^90 36' ^ 120 ' > U 

-(A i_ r 2 _|_ J_ r 4 _ J_ 6N r7 

Vl40 240 ~ 72' 720 

T rM ^ 2 + 1 r ^)S 8 + 0(S 9 ) (8) 

T'l 560 1440 ^ 480 5040 > u ~ ^ V" J ' v ^ 

Numerical eigenanalysis of the diffusion equation (J16)) reported in Section 0] 
confirms the high order accuracy and stability of the resultant integration 
scheme with patch boundary conditions from the above operator. 



3 Achieve high order consistency 

Here we demonstrate analytically that appropriate patch boundary condi- 
tions achieve high order consistency for a wide class of pdes. Consider the 
linear pde 

du d 2 u du ^d 3 u d 4 u 

dt dx 2 dx dx 3 dx A ' 
for some constants a, b and c — we have chosen time and space scales so 
that the coefficient of the diffusion term is 1. Specially crafted boundary 
conditions on small patches ensures macroscopic consistency. 
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Following the dynamical systems approach of holistic discretisation [10, 
ITT] we introduce the parameter 7 to control the coupling between patches: 
when 7 = the patches are uncoupled to provide a base for us to apply centre 
manifold theory; but when we subsequently set 7 = 1 we recover a dynamical 
model for the original pde. Modify the operator (JBJ) to invoke the patch 
boundary condition (ptBc) that on x = Xj (noting that the E ±r implies the 
left-hand side is evaluated on the edge of the patch at x = Xj ± rH): 

E ±r Hd xVj = { 1 [^±r5 2 ]+ 1 2 [-{\-\r 2 )^Tr{^~\r 2 )5 i ] 

+ 7 3 [Hi - V + ^ V* ± r(£ - ±r 2 + &r*)S°] 

+ 7 4 Hllo-^ 2 + ^ 4 -^ 6 )M 7 

T^-^ 2 + ^ 4 -^V]}f/ r (10) 

these PtBCs that when 7 = the small patches are decoupled and the re- 
sulting insulating boundary conditions, E ±r d x Vj = , cause the dissipative 
dynamics of Q in each patch to decay exponentially quickly to some con- 
stant field in each patch, namely Vj(x,t) — > Uj for each of the m patches. 1 
For non-zero coupling parameter 7 the subgrid scale patch field is no longer 
constant, and each patch grid value Uj evolves because of the coupling with 
its neighbours. We construct a series solution of the pde © in the coupling 
parameter 7: the first order expression for the microscopic subgrid scale field 
is straightforward, namely 

Vj = U 3 +1 + ±f 2 5 2 ) Uj + C1 H {-\r 2 i + \e) S 2 U 3 + O ( 7 2 , a 2 + b 2 + c 2 ) , 

' (11) 

1 Why are we only using one pair of boundary conditions for the apparently fourth 
order pde One answer is that the PDE may be viewed as an equivalent pde for 

a microscopic simulator that only requires one pair of boundary conditions. For exam- 
ple, if the microscopic simulator is simply a fine scale discretisation of a pde, such as 
u t = —c(5/h)u + (<5 2 /h 2 )u , then only one pair of boundary conditions are needed for the 
simulator, but it has a high order equivalent pde such as ||SJ|. Another answer is that there 
is no physical boundary at the edge of a patch and so we only need resolve smooth subgrid 
fields and for smooth solutions we need only treat the higher order terms as perturbations; 
see that our error terms are expressed as 0(a 2 + b 2 + c 2 ) for this reason. 
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where the microscopic variable £ = (x — Xj)/H ranges over |£| < r ; accom- 
panying these subgrid fields the grid values Uj evolve according to standard 
second order discretisation (upon putting 7 = 1) 



Uj = 7 (J^PUj - |^/) + 0( 7 2 , a 2 + b 2 + c 2 ). 



(12) 



See that the powers of the coupling parameter 7 in the PtBC (jl())l are chosen 
so that discarding terms of 0(^j p ) results in a discrete model, such as (fT2"|) . 
which is of width 2p — 1 in the grid values Uj] for example, the above model 
only involves Uj and Uj±\ . Centre manifold theory e.g.] asserts that for 
small enough 7 all neighbouring solutions are exponentially quickly attracted 
to the resultant model which faithfully describes the dynamics of the system. 
Although no proof is yet available, we anticipate that the case of interest, 
when 7 = 1 , is small enough for this novel theoretical support to still hold. 

In the interim we demonstrate high order consistency. We obtain models 
that resolve more detail of the subgrid microscopic dynamics and its inter- 
action with neighbouring patches by determining higher order terms in the 
coupling parameter 7. Iteration j^j straightforwardly generates higher order 
approximations. 2 For example, discarding terms 0(7 3 ) the subgrid field in 
each patch is modified from (jll|) to 



1 + 7 + ies 2 } + y m 3 - o^ 3 + ue - ^ 4 ] 
+ ch [( 7 - 7 2 )K 3 5 2 + - u^] + 



+ r 2 



H 

-cff(7 - 7 2 )!^ 2 - cH^e - i£)5 4 - -^7 2 ^ 4 



+ r *cH\^U 3 + 0( 1 \a 2 + b 2 + c 2 ). 



(13) 



2 We use the reduce computer algebra package which has free demonstration ver- 
sions available via http://reduce-algebra.com. Obtain our code for this problem from 
http : //www. sci .usq. edu. au/staf f /aroberts/linpbc . red 
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The first line in (113)1 contains the leading few terms in a universal subgrid 
structure for symmetric operators. However, odd operators, such as the 
advection and the dispersion , generate nontrivial subgrid structures 
in each patch, such as those in the second line of (|13|) . which reflect subgrid 
scale interaction of processes. The third and fourth line of the approximate 
field (TH?|) depend upon the patch size r = h/{2H). But physically the 
subgrid scale field in each patch should be independent of the patch size r. 
Although there is some dependence in these approximations, higher orders 
in the coupling parameter 7 remove it. For example, at the beginning of 
the third line in (fT3J) see the term — cif (7 — 7 2 )|£5 2 disappears when we 
set 7 = 1 for the physically relevant approximation. Similarly, computing 
the next order terms in coupling parameter 7 generates terms, in 7 s , which 
cancel the r dependent terms in the third and fourth line of the subgrid 
field (fT3|) . Thus higher order models push any undesirable r dependence to 
higher orders, thereby usefully predicting a subgrid field largely independent 
of the patch size r. 

Simultaneously with the derivation of the subgrid field (j!3|) we determine 
the corresponding evolution of the macroscopic grid values Uj for the pde Q . 
Computing to higher order in the coupling parameter 7 produces refinements 
to the basic discretisation (JT2J); for example, here we discard terms C(7 4 ) to 
determine 

- (tV - hV B ) U, - £ ( 7 V - h 3 * 6 ) u 3 

+ 0{ 1 \a 2 + b 2 + c 2 ). (14) 

Set 7 = 1 to recover a model for the pde (JHJ) supported by centre manifold 
theory. Note how truncating the expansion to different powers of coupling 
parameter 7 changes the width in Uj of the discrete model. With the patch 
boundary conditions (|10j) the model is independent of the patch size r. 

As well as the novel dynamical systems support of exponentially quick 
attractiveness and long term fidelity at finite grid size H, as mentioned ear- 
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lier, another way to assess the model's relevance is to compare the original 
pde with the equivalent pde obtained from model (|14|) in the limit as the 
macroscopic spacing H — > . From (J 14)) , straightforward algebra 3 deduces 
the equivalent pde 



dU 
~dt ~ 

+ H 2 
+ H 4 



d 2 U ®L_ 2 b d3 u 

^ dx 2 ^ dx ^ dx 3 

1 d 4 U i d 3 U 



7 2 a 



d 4 U 



(7 - 7 2 ) 



12 



dx 4 



6 °dx 3 



dx A 
- (7 2 - 7 3 ) 



4 dx 5 



la 



d 6 U 
dx 6 



60 7 72 7 + 90 7 ) y dx6 <*c dx5 j 

+ 0(H 6 ,^,a 2 + b 2 + c 2 ). (15) 



'i~2 _ A 3 
^80 ' 24 ' 



When the coupling parameter 7 = 1 the second and third lines in the equiv- 
alent pde (fTHjl disappear and consequently the diffusion and advection is 
modelled with errors of O (H 6 ), whereas the dispersion and the fourth-order 
dissipation is modelled with errors O (-ff 4 ) • Should you truncate the discreti- 
sation (|14|) to lower orders in coupling parameter 7, there is less cancellation 
in the equivalent pde and the errors are consequently larger. Conversely, the 
errors move to progressively higher orders as more terms in the coupling pa- 
rameter 7 are retained in the centre manifold discretisation 1)14)1 . Our patch 
boundary conditions fllOl) seem to create excellent discretisations for pdes. 



4 The diffusive model is numerically stable 

Although the PtBCs (|Tn|) form consistent models we need to confirm they 
are numerically stable. Indeed many other forms of PtBCs were tried before 
finding one that was both consistent and numerically stable. In this section 

3 See our reduce code from the internet. 
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Table 1: Growth rates A of perturbations from steady state u = 0: for 
diffusion (J16|) with m patches; with gap to patch ratio r — 0.1; n — 11 
points in the microscale grid; and with the fourth order PtBC fll7j1 . 



m 


1 


2,3 


4,5 


6,7 


m + 1 : 2m 


4 


2 x 1(T 12 


-0.946817 


-2.170942 


n/a 


-99.79 


8 


5 x 1(T 12 


-0.996139 


-3.787268 


-7.132829 


-399.1 


16 


2 x 1(T 10 


-0.999758 


-3.984556 


-8.834269 


-1596. 


32 


-2 x 1(T 10 


-0.999987 


-3.999031 


-8.988851 


-6386. 



we explore the gap-tooth simulations of the simple diffusion equation 

du d 2 u 

— = 7—^7 , and 27r-periodic in x. (16) 
ot ox 2 

Imagine we only have access to the dynamics through a microscopic simulator 
of the diffusion (|16|) . here coded by a fine discretisation on n grid points in 
a patch of microscopic size h = rH and with some microscopic time step, 
typically At = lO^-lO" 4 . 

Firstly we implement the PtBC that on the edge of each patch the fine 
discretisation has boundary condition 

[fi5 ±r5 2 -(l- \r 2 )^ =F — \r 2 )5 A ] Uj = Hd x Vj at x = X, ± rH . 

(17) 

Obtain this from the first few terms of (JBJ) or equivalently from PtBC (fTT1|) by 
discarding 0(y 3 ) terms. For the j'th patch this PtBC involves macroscopic 
grid values L^_ 2 , . . . , Uj +2 only. Then systematically perturbing each and 
every microscopic value from zero, there are mn such microscopic values, we 
numerically determined the map of one microscopic time step. 4 Transform 
the eigenvalues \i of this map to growth rates A = log(/x) / At . The mn growth 
rates fall into n groups of m modes. Each group corresponds to a microscopic 

4 In general, the dominant eigenvalues of the time-stepper map may be obtained via 
a matrix-free Krylov subspace iteration Thus for particle simulations we do not 

necessarily need access to all the fine details of the microscale. 
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Table 2: Growth rates of perturbations from steady state u = as for Tabled 
but fewer points in the fine grid, namely n = 7 . 



m 


1 




2,3 


4,5 


6,7 


m + 1 : 2m 


4 


8 x 


io- 13 


-0.947206 


-2.173003 


n/a 


-99.30 


8 


-8 


x 10- 12 


-0.996246 


-3.788826 


-7.138379 


-397.2 


16 


-1 


x KT 11 


-0.999785 


-3.984985 


-8.836383 


-1588. 


32 


8 x 


io- 11 


-0.999994 


-3.999139 


-8.989397 


-6355. 



internal mode of the dynamics, roughly exp(A^t) cos[£n(x — Xj + h/2)/h] 
for growth rate A^ ~ —£ 2 Tc 2 /h 2 for £ = 0, 1, . . . , n — 1 . For £ > 1 these 
are the rapidly decaying microscopic modes internal to each patch seen in 
the initial instants of the simulations of Figures Q and [21 The other group 
of m modes, £ = , with small growth rates, correspond to the relatively 
slowly evolving macroscopic modes of interest that arise through the coupling 
between patches of the microscopic dynamics. Table Q shows the leading 
seven growth rates, and the magnitude of the £ = 1 internal growth rates, 
for various numbers of patches, m = 4,8, 16, 32 . The exact growth rates 
of the diffusion pde ()16|) are A = — k 2 for integer k. See in the table that 
as the number of patches double, the accuracy of the growth rates of the 
macroscopic modes improves by a factor of about 16. This is consistent with 
an (9 (if 4 ) method as predicted for diffusion with PtBC (I17j) . 

Second, we repeat the analysis for fewer subgrid points so that the micro- 
scopic dynamics are not resolved as well. Tabled shows the leading eigenval- 
ues for n = 7 points in each patch. There is no significant difference between 
Tables Hand El indicating that the microscopic resolution, the only difference 
between the two, has little impact on the macroscopic dynamics. No growth 
rate is significantly positive showing the numerical method is stable — the 
leading growth rate is close to zero corresponding to conservation of mate- 
rial. The other dominant growth rates rapidly approach those for diffusion. 

Lastly, consider the diffusive dynamics of (fT7)|) when connected by the 
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Table 3: Growth rates of perturbations from steady state u = as for Tabled 
but with the sixth order PtBC (jTHjl . 



m 


1 


2,3 


4,5 


6,7 


m + 1 : 2m 


4 


8 x 10~ 12 


-0.982238 


-2.457648 


n/a 


-99.79 


8 


4 x KT 11 


-0.999677 


-3.928952 


-7.843254 


-399.1 


16 


4 x 10~ n 


-1.000006 


-3.998708 


-8.967122 


-1596. 


32 


-2 x 1(T 10 


-1.000003 


-4.000023 


-8.999625 


-6386. 



PtBC that at x — Xj ± rH 

Hd x Vj = [^±r5 2 -{\-\r 2 )^Tr{±-\r 2 )5 i 

+ (£ - K + ^ 4 )/^ 5 ± - sir 2 + ik 4 )* 6 ] ^ • (is) 

Obtain this PtBC from the first six terms of (JBJ) or equivalently from PtBC ()10|) 
by discarding 0(7 4 ) terms. Table El demonstrates that the resultant numer- 
ical scheme is stable and has sixth order consistency for the diffusion equa- 
tion. Further, it is these PtBCs we used to simulate the nonlinear dynamics 
of Burgers' equation ((TJ) to create Figures [T] and El 

5 Conclusion 

We achieve higher order accuracy in the gap-tooth scheme using carefully 
crafted patch boundary conditions ( PtBCs). Analytic approximations and 
analysis of numerical steps in time confirm the PtBCs (|17|) and (|18|) are effec- 
tive. Importantly the PtBCs (JT7J) and (JTHJ) do not depend upon the particular 
pde being simulated, thus the PtBCs should work effectively for particle sim- 
ulations for which we do not have an algebraic microscale closure. 

Further, although the predicted microscopic subgrid scale fields do have 
some dependence upon the patch size r, the dependence weakens, by being 
pushed to higher orders in r, in using higher order accuracy patch boundary 
conditions. 
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As shown in Figures Q and El the PtBCs we recommend here appear to 
work well even for the nonlinear dynamics of Burgers equation (JTJ. 
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